%% IMPORT and Organize DATA

names = fullnames_data(choose_variab)
bettername = betternames_data(choose_variab) 
region = region_data(choose_variab)
names_exog = fullnames_data(choose_exog)
Y = data_all(:,choose_variab);
Y_exog=data_all(:,choose_exog);

% chose names for IRFs
names_IRFs = fullnames_data(choose_variab);
transformIRF_used = transformIRFs_data(choose_variab);
labelIRF_used = labelIRFs_data(choose_variab);
Numbers = 1:length(day);
Numbers=Numbers';

% reduce to period for which all measures are available
cond = isnan(Y);
cond2 = sum(cond,2) == 0;
Y = Y(cond2,:);
Y_exog = Y_exog(cond2,:);
year_available     = year(cond2);
month_available    = month(cond2);
day_available    = day(cond2);
Numbers_available = Numbers(cond2);
% plot(Y(:,2));

data_all_available = data_all(cond2,:);
% keep only after a certain period, in case data are available for all variables added for before that period

% full sample
from = 1 
until = 30000; 

cond = and(Numbers_available >= from,Numbers_available <= until);
sum(cond)
Y = Y(cond,:);
Y_exog = Y_exog(cond,:);

year_available = year_available(cond);
month_available = month_available(cond);
day_available    = day_available(cond);
Numbers_available = Numbers_available(cond);

disp(['-->> Data available for :    ', num2str(year_available(1)) 'M' num2str(month_available(1)) 'D' num2str(day_available(1)) ...
                            ' until ', num2str(year_available(end)) 'M' num2str(month_available(end)) 'D' num2str(day_available(end))] );
% ensure Y enters horizontally
Y = Y';
Y_exog = Y_exog';
k = size(Y,1)
T = size(Y,2)
assert(  T > k)
assert( length(year_available) == T)      

year_available_help=zeros(T,1);
for i=1:length(year_available)
year_available_help(i,:) = year_available(i,:) + 0.001*i;
end


%% Estimate RVAR
maxlags = 20;
 
constantD   = 1;
trendD      = 0;
trendBreaks = [];
seasonD     = 0;
moreBreaks  = [];


[Ahat, Resid, Sigmahat, Ahat_determ, X_determ, T_used, df_lost] = ...
         zFC_VAR_Estimation_exog(Y,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog);        
        
year_available_used     = year_available(p+1:end);
assert(length(year_available_used) == size(Resid,2))
month_available_used    = month_available(p+1:end);
day_available_used    = day_available(p+1:end);
Numbers_available_used = Numbers_available(p+1:end);

%% IRFs, point estimation
T_step = [1:1:T_irf]; % used to plot IRFs
shock_considered = 1; % variable shocked

Instrument = data_all(:,choose_shock);
cond_i = isnan(Instrument);
cond2_i = sum(cond_i,2) == 0;
Instrument = Instrument(cond2_i,:);
year_available_ins     = year(cond2_i);
month_available_ins    = month(cond2_i);
day_available_ins = day(cond2_i);
Numbers_available_ins = Numbers(cond2_i);

% select subperiod common to shocks used and variables used
from2_year         = max(year_available_used(1),year_available_ins(1));
from2_month          = max(month_available_used(1),month_available_ins(1));
from2_day          = max(day_available_used(1),day_available_ins(1));
from2          = max(Numbers_available_used(1),Numbers_available_ins(1));
% % % % %from2          = timeaxis_used(59);
until2_year         = min(year_available_used(end),year_available_ins(end));
until2_month         = min(month_available_used(end),month_available_ins(end));
until2_day         = min(day_available_used(end),day_available_ins(end));
until2         = min(Numbers_available_used(end),Numbers_available_ins(end));
Numers_available_used2 = [from2:until2];
%                         
disp(['-->> Period effectively used in the identification:    ', num2str(year_available_ins(1)) 'M' num2str(month_available_ins(1)) 'D' num2str(day_available_ins(1)) ...
                            ' until ', num2str(year_available_ins(end)) 'M' num2str(month_available_ins(end)) 'D' num2str(day_available_ins(end))] );
                      
% Take only the instruments in this subperiod
cond1            = Numbers_available_ins >= from2 & Numbers_available_ins <= until2;
Instruments_used = Instrument(cond1,:);
n_instrum=size(Instrument);
n_instrum=n_instrum(1,2);

% Take only the residuals in this subperiod
cond2      = Numbers_available_used >= from2 & Numbers_available_used <= until2;
Resid_used = Resid(:,cond2);
AAA=[Resid_used' Instruments_used];
%assert(length(timeaxis_used2) == size(Resid_used,2))

names_instrument = fullnames_data(choose_shock) ;
[bla,g] = size(Instruments_used) ;

% Identify the model
% choose correct position for indicator variable
position_indicator = 1;


[imp_vect_relat, absolute_impulse_vector, ...
             F_stat_test1, p_values_test1, signif_test1, Rsquared_test1, Betas_test1, ...
             F_stat_test2, p_values_test2, signif_test2, Rsquared_test2, ...
             correl_instrument_shock, a_vect, estim_shocks] = zFC_IdentifExtInstruments5(Resid_used, Instruments_used, position_indicator, Sigmahat);
% 
table_validity_rows = char('beta', 'numb. of stars', 'numb. of obs.', 'F statistic', 'R squared', 'obs. of the proxy used');
table_validity_rows = cellstr(table_validity_rows); 
table_validity = [Betas_test1'; signif_test1(:,1)'; length(Instruments_used)*ones(1,k)/1000; F_stat_test1'; Rsquared_test1']    

nonzero_instrument = nonzeros(Instruments_used);
n_nonzero_instrument = length(nonzero_instrument)


B_extin    = [absolute_impulse_vector, zeros(k,k-1)];

if sd_shock == 1
    shock_size_extin = 1;
else
    shock_size_extin = shock_size/absolute_impulse_vector(1);
end
        
IRF_extin  = zFC_IRFs(Ahat, B_extin, T_irf, 1, shock_size_extin);


%% Bootstrap on the reduced form model 

if bootstrap_on == 1 
    
    Yinit           = Y(:,1:p);
    Ahat_loop       = NaN*ones(k,k*p,Loops);
    Resid_used_loop = NaN*ones(k,length(Resid_used),Loops);
    Sigmahat_loop   = NaN*ones(k,k,Loops);
    Instruments_used_loop = NaN*ones(length(Instruments_used),n_instrum,Loops);

    for l = 1:Loops

        % Generate a vector that flips signs randomly
        step = randn(T_used,1);
        change_sign = -1*ones(T_used,1);
        change_sign(step > 0) = 1;

        % Generate pseudo data
        Resid_pseudo = Resid*diag(change_sign);
        % Resid_pseudo = Resid;      % like this it should replicate exactly the data

        Y_pseudo = zFC_GenerateVARData(Yinit, Ahat, Resid_pseudo, Ahat_determ, X_determ);

        % Estimate the VAR on the pseudo data
        [Ahat_loop(:,:,l), Resid_step, Sigmahat_loop(:,:,l), ~, ~, ~, ~] = ...
                       zFC_VAR_Estimation_exog(Y_pseudo,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog);

        % Store residuals for the subpart that overlaps with the shocks used
        Resid_used_loop(:,:,l) = Resid_step(:,cond2);

        % Flip sign of the instruments using the same order in the bootstrapped residuals, consider only the part that overlaps with the residuals
        Instruments_used_loop(:,:,l) = diag(change_sign(cond2))*Instruments_used;
    end

    % IRFs, error bands
    IRF_chol_loop  = NaN*ones(k,T_irf,Loops);
    IRF_extin_loop = NaN*ones(k,T_irf,Loops);
    IRF_extin_new_loop = NaN*ones(k,T_irf,Loops);

    for l = 1:Loops

        % External instrument, using only monthly data
        [~, ~, relative_impulse_vector_loop, ~, incons_impulse_vector_loop] = zFC_IdentifExtInstruments2(Resid_used_loop(:,:,l),Instruments_used_loop(:,:,l));
        b11_loop = zFC_ScaleInExternalInstrument(Sigmahat_loop(:,:,l), relative_impulse_vector_loop);
        absolute_impulse_vector_loop = relative_impulse_vector_loop*b11_loop;

        B_extin_loop    = [absolute_impulse_vector_loop, zeros(k,k-1)];
        
        if sd_shock == 1
            shock_size_extin_loop = 1; 
        else
         shock_size_extin_loop = shock_size/absolute_impulse_vector_loop(1);
        end

        IRF_extin_loop(:,:,l)  = zFC_IRFs(Ahat_loop(:,:,l), B_extin_loop, T_irf, 1, shock_size_extin_loop);

    end

    % Compute error bands

    IRF_extin_bands = NaN*ones(k,T_irf,2);
    interval       = 10;
    for ttt = 1:T_irf 
        for kkk = 1:k
            if cholesky_on == 1
                IRF_chol_bands(kkk,ttt,1)   =  prctile(IRF_chol_loop(kkk,ttt,:),interval/2);
                IRF_chol_bands(kkk,ttt,2)   =  prctile(IRF_chol_loop(kkk,ttt,:),100-interval/2);   
            end
            IRF_extin_bands(kkk,ttt,1)  =  prctile(IRF_extin_loop(kkk,ttt,:),interval/2);
            IRF_extin_bands(kkk,ttt,2)  =  prctile(IRF_extin_loop(kkk,ttt,:),100-interval/2);   
        end
    end

end
%% Final figure

figure(2) 

for i = 1:k
    subplot(rows_fig,cols_fig,i)

    if bootstrap_on == 1
    
    fill([T_step,fliplr(T_step)],[transformIRF_used(i)*IRF_extin_bands(i,:,2), fliplr(transformIRF_used(i)*IRF_extin_bands(i,:,1))], ...
                        grey,'EdgeColor','none', 'FaceAlpha', 0.5), hold on
    end
    
     a = plot(T_step, transformIRF_used(i)*IRF_extin(i,:),     'k-', 'Linewidth', 2);  hold on

     if cholesky_on == 1
         fill([T_step,fliplr(T_step)],[transformIRF_used(i)*IRF_chol_bands(i,:,2), fliplr(transformIRF_used(i)*IRF_chol_bands(i,:,1))], ...
                            blue,'EdgeColor','none', 'FaceAlpha', 0.5), hold on
          b = plot(T_step, transformIRF_used(i)*IRF_chol(i,:),     'b--', 'Linewidth', 2), hold on
     end
    title_i = strcat(bettername(i)) ;
    plot(T_step,zeros(T_irf,1), 'k')
    title(title_i,'FontSize',font_size)
       ylabel(labelIRF_used(i),'FontSize',font_size)
    axis tight
    set(gca,'box','off','FontSize',font_size)
    set(gca,'box','off')
end

%% FOR VD AND HD
% Estimate the a vector and the structural shocks 
Sigma11 = Sigmahat(1,1);
Sigma21 = Sigmahat(2:end,1);
Sigma22 = Sigmahat(2:end,2:end);
b11 = absolute_impulse_vector(1);
b21 = absolute_impulse_vector(2:end);
b12primeB22inv = (Sigma21 - b11*b21)'*(Sigma22 - b21*b21')^(-1);
a_vect_prime = (b11 - b12primeB22inv*b21)^(-1)*[1, -b12primeB22inv];
a_vect_final = a_vect_prime';
estim_shocks_final = a_vect'*Resid;

estim_shocks_final_cum = NaN*ones ;
for t = 1:T_used
    estim_shocks_final_cum(t) = sum(estim_shocks_final(1:t));
end

Y_used = Y(:,p+1:end);
Autoregressive = Ahat;
B = [absolute_impulse_vector, randn(k,k-1)];
StructShocks = [estim_shocks_final; randn(k-1,T_used)];


%% VARIANCE DECOMPOSITION

if fevd_on == 1 
    
    horizons = [1,5,10,50,100]';
    H = length(horizons);

    Var_struct_shocks = eye(k); % normalization used both in Choleski and in external instruments

    % compute moving average representation
    [MovAverage_reduced, ~] = zFC_RewriteVAR(Ahat, randn(k,k), max(horizons));

    B_extin2 = [absolute_impulse_vector, randn(k,k-1)];
    [~, MovAverage_structural_extin] = zFC_RewriteVAR(Ahat, B_extin2, max(horizons));

    Vdecomp  = NaN*ones(H,k);
    for h = 1:H
        hor = horizons(h);    
        for variab = 1:k
            shock = 1;
            Vdecomp(h,variab) = zFC_VarianceDecomposition_partially(shock, variab, hor, Sigmahat, MovAverage_reduced, Var_struct_shocks, MovAverage_structural_extin);   
        end
    end

    table_FEVD = [horizons, round(Vdecomp,1)]

end

%% HISTORICAL DECOMPOSITION

if hd_on == 1 
    
    t_backuntil = 1; 
    t_until = Numbers_available_used(end)-Numbers_available_used(1)+1;
    tick=Numbers_available_used;

    [Y_datavalues, Base_all, Y_fromeachshocks_all] = zFC_HistoricalDecomposition_partially(Y_used, Autoregressive, B, Resid, StructShocks, t_backuntil, t_until);

    figure(6)

    subplot(2,1,1)
       plot(Numbers_available_used,exp(Y_datavalues(2,:)), 'k', 'Linewidth', 2), hold on
       plot(Numbers_available_used,exp(Y_datavalues(2,:)' - squeeze(Y_fromeachshocks_all(i,1,:))), '-- r', 'Linewidth', 1.5), hold on
       axis tight
       title('Historical decomposition','FontSize',12)   
       set(gca,'XTick',[min(tick),min(tick)+2000,min(tick)+4000,max(tick)])
       set(gca,'XTickLabel',[year(min(tick),:),year(min(tick)+2000,:),...
       year(min(tick)+4000,:),year(max(tick),:)])
       set(gca,'box','off','FontSize',10)
       h=legend('USD/JPY', 'USD/JPY without Interv.') ; 
       set(h,'Location','northwest')

    subplot(2,1,2)
        plot(tick, estim_shocks_final,'k','Linewidth', 1), hold on
        plot(tick, estim_shocks_final_cum,'k','Linewidth', 2), hold on
        title('Intervention shocks','FontSize',12)   
        h= legend('Estimated shocks', 'Cumulated shocks') ;
        set(h,'Location','northwest')
        axis tight
        set(gca,'XTick',[min(tick),min(tick)+2000,min(tick)+4000,max(tick)])
        set(gca,'XTickLabel',[year(min(tick),:),year(min(tick)+2000,:),...
        year(min(tick)+4000,:),year(max(tick),:)])

end
